Fast, adaptive, high order accurate discretization of the Lippmann-Schwinger equation in two dimension

We present a fast direct solver for two dimensional scattering problems, where an incident wave impinges on a penetrable medium with compact support. We represent the scattered field using a volume potential whose kernel is the outgoing Green's function for the exterior domain. Inserting this representation into the governing partial differential equation, we obtain an integral equation of the Lippmann-Schwinger type. The principal contribution here is the development of an automatically adaptive, high-order accurate discretization based on a quad tree data structure which provides rapid access to arbitrary elements of the discretized system matrix. This permits the straightforward application of state-of-the-art algorithms for constructing compressed versions of the solution operator. These solvers typically require $O(N^{3/2})$ work, where $N$ denotes the number of degrees of freedom. We demonstrate the performance of the method for a variety of problems in both the low and high frequency regimes.Comment: 18 page

An accurate, fast, mathematically robust, universal, non-iterative algorithm for computing multi-component diffusion velocities

Using accurate multi-component diffusion treatment in numerical combustion studies remains formidable due to the computational cost associated with solving for diffusion velocities. To obtain the diffusion velocities, for low density gases, one needs to solve the Stefan-Maxwell equations along with the zero diffusion flux criteria, which scales as $\mathcal{O}(N^3)$, when solved exactly. In this article, we propose an accurate, fast, direct and robust algorithm to compute multi-component diffusion velocities. To our knowledge, this is the first provably accurate algorithm (the solution can be obtained up to an arbitrary degree of precision) scaling at a computational complexity of $\mathcal{O}(N)$ in finite precision. The key idea involves leveraging the fact that the matrix of the reciprocal of the binary diffusivities, $V$, is low rank, with its rank being independent of the number of species involved. The low rank representation of matrix $V$ is computed in a fast manner at a computational complexity of $\mathcal{O}(N)$ and the Sherman-Morrison-Woodbury formula is used to solve for the diffusion velocities at a computational complexity of $\mathcal{O}(N)$. Rigorous proofs and numerical benchmarks illustrate the low rank property of the matrix $V$ and scaling of the algorithm.Comment: 16 pages, 7 figures, 1 table, 1 algorith

New algebraic fast algorithms for NN-body problems in two and three dimensions

This article presents two new algebraic algorithms to perform fast matrix-vector product for $N$-body problems in $d$ dimensions, namely nHODLR$d$D (nested algorithm) and s-nHODLR$d$D (semi-nested or partially nested algorithm). The nHODLR$d$D and s-nHODLR$d$D algorithms are the nested and semi-nested version of our previously proposed fast algorithm, the hierarchically off-diagonal low-rank matrix in $d$ dimensions (HODLR$d$D), respectively, where the admissible clusters are the certain far-field and the vertex-sharing clusters. We rely on algebraic low-rank approximation techniques (ACA and NCA) and develop both algorithms in a black-box (kernel-independent) fashion. The initialization time of the proposed hierarchical structures scales quasi-linearly. Using the nHODLR$d$D and s-nHODLR$d$D hierarchical structures, one can perform the multiplication of a dense matrix (arising out of $N$-body problems) with a vector that scales as $\mathcal{O}(pN)$ and $\mathcal{O}(pN \log(N))$, respectively, where $p$ grows at most poly logarithmically with $N$. The numerical results in $2$D and $3$D $(d=2,3)$ show that the proposed nHODLR$d$D algorithm is competitive to the algebraic Fast Multipole Method in $d$ dimensions with respect to the matrix-vector product time and space complexity. The C++ implementation with OpenMP parallelization of the proposed algorithms is available at \url{https://github.com/riteshkhan/nHODLRdD/}.Comment: 32 page

Fast and scalable Gaussian process modeling with applications to astronomical time series

The growing field of large-scale time domain astronomy requires methods for probabilistic data analysis that are computationally tractable, even with large datasets. Gaussian Processes are a popular class of models used for this purpose but, since the computational cost scales, in general, as the cube of the number of data points, their application has been limited to small datasets. In this paper, we present a novel method for Gaussian Process modeling in one-dimension where the computational requirements scale linearly with the size of the dataset. We demonstrate the method by applying it to simulated and real astronomical time series datasets. These demonstrations are examples of probabilistic inference of stellar rotation periods, asteroseismic oscillation spectra, and transiting planet parameters. The method exploits structure in the problem when the covariance function is expressed as a mixture of complex exponentials, without requiring evenly spaced observations or uniform noise. This form of covariance arises naturally when the process is a mixture of stochastically-driven damped harmonic oscillators -- providing a physical motivation for and interpretation of this choice -- but we also demonstrate that it can be a useful effective model in some other cases. We present a mathematical description of the method and compare it to existing scalable Gaussian Process methods. The method is fast and interpretable, with a range of potential applications within astronomical data analysis and beyond. We provide well-tested and documented open-source implementations of this method in C++, Python, and Julia.Comment: Updated in response to referee. Submitted to the AAS Journals. Comments (still) welcome. Code available: https://github.com/dfm/celerit

HODLRddD: A new Black-box fast algorithm for NN-body problems in dd-dimensions with guaranteed error bounds

In this article, we prove new theorems bounding the rank of different sub-matrices arising from these kernel functions. Bounds like these are often useful for analyzing the complexity of various hierarchical matrix algorithms. We also plot the numerical rank growth of different sub-matrices arising out of various kernel functions in $1$D, $2$D, $3$D and $4$D, which, not surprisingly, agrees with the proposed theorems. Another significant contribution of this article is that, using the obtained rank bounds, we also propose a way to extend the notion of \textbf{\emph{weak-admissibility}} for hierarchical matrices in higher dimensions. Based on this proposed \textbf{\emph{weak-admissibility}} condition, we develop a black-box (kernel-independent) fast algorithm for $N$-body problems, hierarchically off-diagonal low-rank matrix in $d$ dimensions (HODLR$d$D), which can perform matrix-vector products with $\mathcal{O}(pN \log (N))$ complexity in any dimension $d$, where $p$ doesn't grow with any power of $N$. More precisely, our theorems guarantee that $p \in \mathcal{O} (\log (N) \log^d (\log (N)))$, which implies our HODLR$d$D algorithm scales almost linearly. The $\texttt{C++}$ implementation with \texttt{OpenMP} parallelization of the HODLR$d$D is available at \url{https://github.com/SAFRAN-LAB/HODLRdD}. We also discuss the scalability of the HODLR$d$D algorithm and showcase the applicability by solving an integral equation in $4$ dimensions and accelerating the training phase of the support vector machines (SVM) for the data sets with four and five features.Comment: 35 pages, 23 figures, 14 table